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We study electronic properties of graphene with finite concentration of vacancies or other resonant 
scatterers by a straightforward lattice Quantum Monte Carlo calculations. Taking into account real¬ 
istic long-range Coulomb interaction we calculate distribution of spin density associated to midgap 
states and demonstrate antiferromagnetic ordering. Energy gaps are open due to the interaction 
effects, both in the bare graphene spectrum and in the vacancy/impurity bands. In the case of 5 % 
concentration of resonant scatterers the latter gap is estimated as 0.7 eV and 1.1 eV for graphene 
on boron nitride and freely suspended graphene, respectively. 
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Defects effect enormously on electronic properties of 
graphene and other Dirac materials. In particular, va¬ 
cancies in graphene are known to create mid-gap states 
[1, 2] which, together with electron-electron interaction, 
can result in appearance of magnetic moments and rich 
many-body phenomena (see the review of early works in 
Refs. [2, 3] and recent experimental and theoretical pa¬ 
pers [4-7]). Note that hydrogen adatoms and some uni¬ 
valent organic admolecules (resonant scattering centers) 
produce a very similar electronic structure [8]. Qualita¬ 
tively, it is explained by the fact that sp 3 state of carbon 
atom originated from its bond with a univalent adatom 
or admolecule makes it unavailable for p z electrons (tt- 
orbitals) at energies at the neutrality point plus minus 
several electronvolts; for these electrons such atom is 
just cut from the lattice. Thus, discussing “vacancies” 
we will keep in mind also these cases; moreover, they 
are even closer to a simple model of vacancy just as a 
missed site in the honeycomb lattice (the model which 
will be used in our calculations) since the real vacancy 
produces very strong lattice distortions effecting essen¬ 
tially on the electronic structure [9]. The case of finite 
concentration of vacancies is quite complicated even at a 
single-particle level [11—13]. Here we consider this case 
taking into account a realistic model of Coulomb interac¬ 
tion in graphene [14] via the straightforward lattice quan¬ 
tum Monte Carlo (QMC) simulations. Keeping in mind 
what was said above, the best and easiest experimen¬ 
tal realization would be partially hydrogenated graphene. 
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We will study the antiferromagnetic (AFM) phase tran¬ 
sition driven by the presence of adatoms. It manifests 
itself in emergent magnetic moments concentrated near 
the adatoms as well as in the band gap opening. Both fea¬ 
tures of this phase transition are very important. Band 
gap, which is controllable via hydrogenation, offers an 
interesting possibility for prospective graphene applica¬ 
tions in semiconductor devices. Emergent magnetic mo¬ 
ments could be one of the reasons for short spin relax¬ 
ation time in graphene, which is an essential obstacle for 
producing efficient spintronic devices [15, 16]. We em¬ 
phasize that unlike the previous density functional theory 
study of vacancies in graphene [9, 10], here we treat very 
large sample with random (non-regular) distribution of 
vacancies using unbiased QMC. Thus we could estimate 
how close vacancies influence each other in various ge¬ 
ometrical configurations (for example, by calculation of 
variations of magnetic moment). Also we could extract 
energy of midgap state associated with any particular 
scatterer and observe variations of these energies from 
one scatterer to another depending on its surroundings. 
These measurements give us an opportunity to estimate 
realistic width of midgap energy band for a fixed dis¬ 
tribution of scatterers. The demonstrated possibility of 
QMC simulations of large samples with arbitrary posi¬ 
tions of scatterers is even more important because in real 
samples adatoms can migrate and form clusters. Now 
this phenomenon is in focus of intensive research [17, 18]. 
Our technique can be applied for these real spatial con¬ 
figurations of scatterers. In principle, one can calculate 
the potential of interaction between adatoms within the 
QMC and to model the formation of clusters which is an 
interesting project for future. Here we restrict ourselves 
only by the case of random distribution of defects. 
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We start from the tight-binding Hamiltonian for non¬ 
interacting electrons with staggered mass term which is 
essential in our simulations for the following reasons: ( 1 ) 
it eliminates zero mode in spectrum of quasiparticles thus 
making the fermionic operator M (see below) invertible; 
( 2 ) it serves as a seed for antiferromagnetic phase transi¬ 
tion which we will study. The initial Hamiltonian with¬ 
out interaction and adatoms reads: 

Htb =-t + h - C -) ± 

<x,y> 

X 

where t = 2.7 eV, the sum goes over all pairs of 

<x,y> 

nearest-neighbor sites of the graphene honeycomb lattice 
(we impose periodic spatial boundary conditions as in 
Refs. [21, 22]) and the mass term has different sign at 
different sublattices. Here a' x a x -[• and a\. j_, a x ,i are 
the creation/annihilation operators for spin up and spin 
down electrons at 7r-orbitals. 

Next, we introduce the electrostatic interaction with 

potentials V xy : He = \ J2Vx V q x q y , where q x = 

x,y 

+ a' x , a x ^ — 1 is the operator of electric charge 
at lattice site x. The whole matrix V xy is constructed 
in the following way: at small distances (on-site interac¬ 
tion ( 14 K = Voo) an d interactions with the nearest (Voi), 
next-to-nearest (V 02 ) and next-to-next-to-nearest (V 03 ) 
neighbours) we use the potentials calculated by the con¬ 
strained RPA method [14]; at larger distances we use 
ordinary Coulomb 14^ = Vcroi/r xy . Parameter Vc de¬ 
fines the strength of the Coulomb tail. Since roi is the 
distance between the nearest neighbours, Vq is equal to 
the Coulomb repulsion energy at the distance of conju¬ 
gated C-C bond. There is a small difference in calibra¬ 
tion of our model (t = 2.7 eV) and the calculations [14], 
where the hopping parameter was 2.8 eV but any possi¬ 
ble changes of observable quantities due to this difference 
are definitely outside the accuracy of the applied QMC 
technique. 

We use two settings of interaction potentials. The first 
one is called “ordinary potentials” and corresponds to 
freely suspended graphene (the dielectric constant of en¬ 
vironment is equal to one). At small distances it is sim¬ 
ply the set of potentials from the paper [14]: Voo = 9.3 
eV, V 01 = 5.5 eV, V 02 = 4.1 eV and V 03 = 3.6 eV. 
The strength of Coulomb tail is defined by the ratio 
V 03 = Vcr 0 i/ro 3 , thus Vc = 7.2 eV. It is a simple contin¬ 
uous extension of the potentials at small distances. The 
second setting is specified as “screened potentials” and 
corresponds to graphene at a substrate with the dielec¬ 
tric constant e = 4.5 which is roughly the value reason¬ 
able for both boron nitride and SiC> 2 . In this case the 
Coulomb tail is (1 + e)/2 = 2.75 times weaker while V 02 
and V 03 are 1.5 times smaller. Voo and V 01 are untouched 
since the screening by substrate should be irrelevant at a 
few interatomic distances. 
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FIG. 1: (a) Antiferromagnetic order parameter at 
various temperatures, calculated on different lattices in 
presence of adatoms and without them, (b) 
Temperature dependence of antiferromagnetic order 
parameter in the zero bare mass limit. Calculation was 
performed on the 12 x 12 lattice. 

As was already mentioned, the vacancy or adatom is 
modeled by setting to zero all hoppings to the vacant site. 
We also exclude the vacant sites from the interaction 
term He because they have constant zero charge. We 
employ Hybrid Monte-Carlo algorithm, broadly used in 
lattice QCD. It was applied earlier for studies of graphene 
in the papers [19, 20], where so-called staggered fermions 
were used to model low-energy effective field theory of 
graphene. This algorithm was developed further in Refs. 
[21, 22, 25]. At first, we perform Suzuki-Trotter decom¬ 
position of partition function exp(— (3H) in order to rep¬ 
resent it as a functional integral over trajectories in Eu¬ 
clidean time. In order to get rid of four-fermionic terms in 
the Hamiltonian, we use Hubbard-Stratonovich transfor¬ 
mation and obtain the following partition function after 
integrating out fermionic fields: 

Tre-^S Jvip x , n e- s ^-^\det(M[tp x , n ])\ 2 , (2) 

where tp x n is the Hubbard-Stratonovich field for times- 
lice n and spatial coordinate x. 5rN t = j3, where St is the 
step in Euclidean time, N t is number of steps and /? is in- 
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FIG. 2: Distribution of average spin. Color scale 
corresponds to (S z ) at the site in the zero bare mass 
limit. 


verse temperature, M is the fermionic operator (inverse 
fermionic Green’s function at a given configuration of 
auxiliary field). We use its particular form [22], In more 
details (including issues with continuous limit St —> 0) it 
was discussed in Ref. [23]. The particle-hole symmetry 
for graphene at neutrality point makes the integration 
weight in (2) positive due to appearance of the squared 
modulus of the determinant, thus, we have no fermionic 
sign problem [24] . For both sets of inter-electron interac¬ 
tion potentials, the action of the Hubbard-Stratonovich 
field S [</5x,n] is also a positive definite quadratic form. 
Thus we generate configurations of by a Monte- 
Carlo method and calculate physical quantities as av¬ 
erages over the generated configurations. Here we follow 
Refs. [21, 22, 25] and use so-called ^-algorithm. 

We used lattices with spatial sizes 18 x 18, 24 x 24 and 
36 x 36 in order to study finite-size effects. We studied 
lattice with 5 % adatoms (in the most of calculations), 
scattered uniformly in the whole sample. Three differ¬ 
ent temperatures were studied: T=0.5 eV (corresponds 
to N t = 20), T=0.125 eV (N t = 80) and T=0.0625 eV 
(N t = 160). For all temperatures we generated config¬ 
urations with four masses, for example in the case of 
T=0.125 eV we used m = 0.05,0.1,0.15,0.2 eV. Physical 
results are obtained via extrapolation to zero mass. In 
all calculations except energies of midgap states we use 
“ordinary potentials”. 

According to Lieb theorem for the Hubbard model 
[2, 26] the ground state for the case of vacancies equally 
distributed between two sublattices should be spin sin¬ 
glet, and there are no physical reasons to expect that the 
long-range character of Coulomb interactions can change 
this conclusion. Keeping in mind that single vacancy or 
adatom induces magnetic moment one should consider 
opportunity of antiferromagnetic ordering at finite con¬ 
centration (ferromagnetism is impossible). In this case 
the order parameter is the difference in average spin be¬ 
tween sublattices (denoted as A and B in the formula): 


( An ) = (wj E E 

x€A x€B 

A lx vi) )> Na and Nb are the overall number of sites in 
A and B sublattice, respectively. The results are pre¬ 
sented in Fig. la. In case of the highest temperature 
(0.5 eV) the order parameter is equal to zero in the phys¬ 
ical limit of zero bare mass disregarding the presence of 
adatoms. Only at lower temperature (0.125 eV) the order 
parameter acquires nonzero value in presence of adatoms 
and remain almost stable with further decreasing of the 
temperature (0.0625 eV). Fig. lb presents the temper¬ 
ature dependence of AFM order parameter. This calcu¬ 
lation was also performed using one particular random 
distribution of adatoms for each concentration. One can 
clearly see a sharp transition at certain critical temper¬ 
ature (Neel temperature). The results were fitted with 
“step function” /(T) = C(1 — tanh(6(T — a))), where 
parameter a gives us the value of the critical tempera¬ 
ture. Thus we have an estimation of effective antiferro¬ 
magnetic coupling (for 1 % of defects) J ~ 0.1 eV. This 
value is two orders of magnitude larger than the one es¬ 
timated from recent experimental data for vacancies in 
graphene [27]. This is an important point showing that 
probably exchange interaction is very sensitive to real 
electronic srtucture (we mentioned in the introduction 
that for the real vacancies it is very strongly affected by 
atomic reconstruction) and that the use of the simplest 
one-band tight-binding model instead of full-electron cal¬ 
culations can be dangerous for the problems related to 
magnetism. Recent density functional calculations of ex¬ 
change interactions in single-site hydrogenated or fluo- 
rinated graphene [28] predict complicated noncollinear 
magnetic ground states, in a sharp contrast with predic¬ 
tions of Lieb theorem for the single-band model (satu¬ 
rated ferromagnetism). This issue requires further inves¬ 
tigations. 

Spatial distribution of electron spin density is pre¬ 
sented in Fig. 2. It represents the quantity f x = 
{a* x at each lattice site. Since the particle-hole sym¬ 
metry is unbroken, the equality {a\. ^a Xj -|-)+ (&][, ^a Xj j.) = 1 
is satisfied exactly for each lattice site. It means that re¬ 
gions with positive f x have non-compensated spin up and 
negative f x corresponds to non-compensated spin down. 
It is clearly seen that antiferromagnetic order is gener¬ 
ated in the vicinity of adatoms. Moreover, one isolated 
adatom has nonzero average spin (see the first row in 
the table I). This spins tend to be parallel for adatoms 
at one sublattice and antiparallcl for adatoms at differ¬ 
ent sublattices. If adatoms are placed equivalently at 
both sublattices, they generate the same spin excess at 
both sublattices and thus the full spin will be close to 
zero. This means that, indeed, the statement of the Lieb 
theorem [26] remains correct in the case of long-range 
Coulomb interaction. More detailed description of spin- 
spin correlations including explicit evidence of sponta¬ 
neous breaking of SU(2) spin symmetry is presented in 
supplementary material. 

In order to characterize the correlation between 
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Region on the map 

M = 

= 2(5.) 

R3 

0.530 

± 

0.016 

R4 

1.330 

± 

0.026 

R5 

1.542 

± 

0.026 

R6 

2.70 

± 

0.04 

Rl 

3.20 

± 

0.04 


TABLE I: Average magnetic moment of different 
configurations of adatoms. R3 corresponds to one 
isolated adatom, R4 corresponds to 2 adatoms at the 
distance of 2 lattice steps, R5 contains 2 adatoms at the 
distance of 1 lattice step, R1 and R6 contain 4 adatoms 
(denser configuration in case of Rl). 



FIG. 3: Energies of midgap states for two sets of 
inter-electron potentials. Each state correspond to one 
isolated vacancy marked with a number in the figure 2. 

Center and width of the bands are calculated in the 
limit to —> 0. Center is average over energies of all 
states in each band, width is equal to doubled 
dispersion. T=0.125 eV. Real physical situation is 
restored in the limit m —> 0. (inset) Energy gap 
between ’’normal” energy bands. All values correspond 
to the K-point in Brillouine zone. 


adatoms at one sublattice quantitatively, we have mea¬ 
sured full magnetic moment M = 2^b(S z ) of some spa¬ 
tial configurations of vacancies. Results (in units of 
Bohr’s magneton) are summarized in the table I. One can 
observe strong dependence of magnetic moment on the 
geometry of the adatoms configurations. For example, 
two adatoms at the distance of 1 lattice step have mag¬ 
netic moment 1.5 times larger than two isolated adatoms. 

The second set of calculations is devoted to measure¬ 
ments of mass gap in the presence of vacancies/adatoms. 
We study two types of energy bands: “normal” energy 
bands which transfer into Dirac cones in the absence of 


adatoms and midgap states which are concentrated in 
the vicinity of isolated adatoms. In the latter case we 
perform calculations for both sets of potentials to mea¬ 
sure the influence of screening on the energies of midgap 
states. Calculation of the energies is based on the two- 
point Euclidean Green functions contracted with some 
guess for projector to a wavefunction ip(x) of the state 
we are interested in: 

C(t) = Tr . (3) 

x ,y 

At large enough r this correlator is proportional to e~ rE ° 
where Eq is energy of the state under study. In the 
case of “normal” energy band we use lattice exponent 
exp ( ikx ) concentrated at one sublattice with wave vec¬ 
tor k at the K-point of Brillouin zone as a guess for 
the wavefunction. Therefore we are able to estimate the 
lower bond of the energy band and the energy gap be¬ 
tween these bands. In the case of midgap state we guess 
that wavefunction is concentrated in the three nearest 
neighbours of the vacant site. In order to check these 
measurements we perform the same calculation for freely 
suspended graphene without vacancies. In this case the 
gap should be equal to zero in the zero bare mass limit 
[22]. Results for the “normal” energy band are presented 
in the inset on Fig. 3. In presence of vacancies we use 
simple linear fit. Without vacancies the polynomial fit 
4>(m) = C0 + C1TO + C2TO 2 is employed. For the largest lat¬ 
tice (36 x 36) Co is zero within the errorbars, so the fitting 
works well and this lattice is large enough to reproduce 
zero gap in the mbare —• > 0 limit. For smaller lattices one 
can observe non zero cq due to large finite-size effects. 

For midgap states we used wavefunctions concentrated 
near 12 relatively isolated adatoms, marked with black 
numbers in the figure 2. Summarizing all these calcula¬ 
tions (see Fig. 3) we conclude that states concentrated 
near adatoms form two rather broad bands in between 
of “normal” energy bands. Positive and negative ener¬ 
gies of midgap states correspond to adatoms at differ¬ 
ent sublattices. The gap between these bands is calcu¬ 
lated as a distance between two levels with the small¬ 
est absolute values of energies. It can exceed 1 eV for 
suspended graphene, but decreases for graphene at sub¬ 
strate. Width of the bands is a measure of the inter¬ 
play between midgap states concentrated near different 
vacancies. Obviously, if concentration of adatoms tends 
to zero, the energies of the midgap states will be almost 
constant. The same effect is observed here in case of sup¬ 
pressed Coulomb tail: midgap states near isolated vacan¬ 
cies recognize their surrounding more poorly. 

To conclude, electron-electron interactions for finite 
concentration of adatoms lead to antiferromagnetic or¬ 
dering, in a qualitative agreement with Lieb theorem de¬ 
spite its formal inapplicability to systems with long-range 
Coulomb interactions. Even probably more interestingly 
they result in gap opening: a “big gap” of the order of 
several eV at the K point and “smaller gap” (but still 
quite noticeable, about 1 eV for 5% concentration of 
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adatoms) in the mid gap states. The latter prediction 
can be checked by measuring optics for chemically func¬ 
tionalized graphene. One could expect that the effects 
of disorder will smear the gap in the defect band trans¬ 
forming it into a pseudogap. One could hope however 
that the two-peak structure characteristic of the pseu- 
dogapped state can be distinguished from a single-peak 
structure centered new the neutrality point predicted by 
the single-particle picture [29] . Alternatively, the recon¬ 
struction of the defect band can be studied experimen¬ 
tally by scanning tunneling spectroscopy. 
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Supplemental material 

More close view at average spin is presented in the 
figures SI and S2. They demonstrate dependence of 
(S 3 2 ) and (<Si 2 ) on bare mass in presence of vacancies 
(S2) and without them (SI). S is a full spin for some 
region of the lattice: Si = YIxgrSx ,», where S x> i = 

l/2(a^ t , Si,;)' 7 * (a*’l) 5 °* are Pauli matrices. sj a) 

( ') 2 

is a spin for particular sublattice a = 1,2. Since (£. ) 

are even in bare mass, we use the following polynomial 
extrapolation: (f>(m) = Co + cim 2 + C 2 m 4 . It is clearly 
seen that both sublattice and SU(2) spin symmetry are 
restored in the limit of zero mass in absence of vacancies 
((S[ a ^ ) = ( S 3 “) ) at both sublattices in the limit m —» 0, 

see SI). Another situation can be observed in presence of 

/ ,2 

adatoms at S2. We show (S) a ) for the lattice with 5 % 
adatoms choosing the region where adatoms are concen¬ 
trated mostly at one sublattice (region Rl in the figure 
2 in the main text). It is clear that both sublattice and 

SU(2) spin symmetry are broken: (S^ ) ^ ( S^ ) at 
both sublattices even in the zero mass limit. 

Now let us turn to spin-spin correlations for different 

spatial configurations of adatoms. Figures S3 and S4 

, , 2 

demonstrate the dependence of (.S',- ) on the size of the 

system. In the region Rl (see figure 2) adatoms are con¬ 
centrated at the 2d sublattice, while in the region R2 
adatoms are placed equivalently at different sublattices. 
In the first case (S3) S^ spin components at different 
sites within one sublattice are correlated with each other: 
(S 3 ’ ; ) ~ N v with v > 1. The strongest correlation 

is within the 1 st sublattice (the “red” one in the fig¬ 
ure 2) while S^ are almost uncorrelated. Again, S[ a ^ 
and Sg Q ^ spin components behave differently because of 
spontaneously broken SU(2) symmetry. Since there is no 



FIG. SI: Dependence of (S 2 ) on bare mass in absence 
of adatoms. Lattice size is 36 x 36, T=0.125 eV. S is 
full spin of a 6 6 cells region of the lattice. Data for the 
second sublattice is not shown since it is exactly 
identical to the first one. 


difference between and S^ a \ the SU(2) symmetry 
is broken up to U(l) rotations. From the calculation of 
full spin (it includes both sublattices) we see that both 
S 3 and Si are anticorrelated at different sublattices, be¬ 
cause v sufficiently decreases in both cases. In the second 
configuration of adatoms (figure S4) S-f - 1 and S are 
correlated equivalently inside one sublattice. All com¬ 
ponents of electron spin at different sublattices are again 
anticorrelated. Thus we really have a spontaneous break¬ 
ing of SU(2) spin symmetry and sublattice Z 2 symmetry 
leading to antiferromagnetic spin ordering. 



Bare mass, eV 


FIG. S2: Dependence of (S 2 ) on bare mass in presence 
of 5% adatoms. Lattice size is 36 x 36, T=0.125 eV. S 
is full spin of the region Rl (see fig. 2 in the main text). 



In N 

FIG. S3: Dependence of (Si 2 ) on the size of the system. 
Calculation is performed on the lattice with 5 % 
adatoms inside the regions Rl (see fig. 2 in the main 
text). T=0.125 eV. All the results are shown in the 
m —> 0 limit. 
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FIG. S4: Dependence of (Si 2 ) on the size of the system. 
Calculation is performed on the lattice with 5 % 
adatoms inside the regions R2 (see fig. 2 in the main 
text). T=0.125 eV. All the results are shown in the 
m —> 0 limit. 





